Asymptotic expansions of the error of

We present and analyse a family of fully discrete spline Galerkin methods for the solution of boundary integral equations with logarithmic kernels. Following the analysis of Galerkin methods of a previous paper, we show the existence of asymptotic expansions of the error. In this work we deal with the numerical solution of the integral equation Z 1 0 log jx( ) x(s)j2g(s)ds = u0(s); 8s; where u0 is a given 1 periodic function and x is a smooth 1 periodic parametrization of the boundary of a simply connected bounded open set in the plane. We will suppose that the logarithmic capacity of that curve is not equal to one, in order to avoid problems of existence and uniqueness of solution. In Crouzeix and Sayas (1995) the existence of asymptotic expansions of the error of Galerkin methods with smoothest splines was shown. However, the practical impossibility of evaluating exactly in most cases the coe cients and the right-hand of the linear system obtained by the Galerkin method demands the use of approximate numerical integration in the process. We will prove that the formulae transforming the Galerkin method into the so-called Galerkin collocation method (see Hsiao, Kopp and Wendland (1980 and 1985)) lead to a formulation of the problem where there exists still an asymptotic expansion of the error in powers of the discretization parameter h. In fact, as in Crouzeix and Sayas (1995), we will be able to prove expansions like kg h Qhg M X k=f(m)hkQhfkk ChM+1; (1) where Qh denotes an interpolation operator onto the discrete space, f(m) is an integer depending on the degree of the splines used as test-trial functions for the numerical method, and g h is the numerical approximation of g. Throughout what follows, C will denote a certain constant, independent of the parameter h, possibly di erent in each occurrence. The symbol O(hk) denotes a function of the parameter h such that divided by hk remains bounded as h goes to zero. 1 The Galerkin method Let us consider a simply connected bounded open set inR2, whose boundary is given by an in nitely often di erentiable 1 periodic non-singular parametrization x. We suppose that the logarithmic capacity of (see Sloan (1992) for example) is not equal to 1. We consider the periodic Sobolev spaces H t equipped with the norm kukt := 0@jû(0)j2 +X k 6=0 jkj2tjû(k)j21A 1 2 ; where û(k) := Z 1 0 u(s)e 2k {sds: We denote by ( : ; : )0 both the inner product in H0, that is, in L2(0; 1), and its extension as the duality bracket between H t and H t for any t. Let V : H t ! H t+1 be the single layer operator, de ned as the continuous extension of V g(s) := Z 1 0 log jx(s) x( )j2g( )d 1 and let b denote the bounded symmetric bilinear form on H 12 b(g; r) := (V g; r)0: To solve approximately the variational problem (P )8><>: nd g 2 H 12 such that b(g; r) = (u0; r)0; 8r 2 H 1 2 ; we consider the space Vh of the periodic smoothest splines of degree m with nodes on the points si+ 12 := i+ 12 h; i 2 Z; where h := 1=N and N is a positive integer. Denoting Ii := (si 12 ; si+ 1 2 ), if m = 0, Vh := fqh 2 H0 : qhjIi 2 P0;8ig; and if m 1 Vh := fqh 2 Cm 1 : qhjIi 2 Pm;8ig; being Ck the space of real 1 periodic k times di erentiable functions and Pk the set of the polynomials of degree lesser than or equal to k. The Galerkin method is the sequence of nite-dimensional variational problems (Ph)8><>: nd gh 2 Vh such that b(gh; rh) = (u0; rh)0; 8rh 2 Vh; uniquely solvable for h small enough. With the application of numerical integration in coe cients and right-hand side of a linear system equivalent to (Ph), this will be transformed into a new problem (P h )8><>: nd g h 2 Vh such that bh(g h; rh) = (u0; rh)h; 8rh 2 Vh; where bh( : ; : ) and ( : ; : )h are respective approximations of the bilinear form b and of the inner product in L2(0; 1), which will be speci ed in the next sections. Let si := hi for all i 2 Z. The interpolation operator Qh maps u 2 C0 to the unique Qhu 2 Vh such that Qhu(zi) = u(zi); 8i; where zi := ( si; if m is even, si 1 2 ; if m is odd. Consequently Qh interpolates on the midpoints of the grid for m even and on the nodes for m odd. 2 Once we have an expansion of bh(Qhg g h; rh) similar to that of Proposition 1 in Crouzeix and Sayas (1995), the uniform inf-sup condition will be used to prove expansions like (1). Since for all rh 2 Vh bh(Qhg g h; rh) = = b(Qhg gh; rh) + (u0; rh)0 (u0; rh)h + bh(Qhg; rh) b(Qhg; rh) (2) the study of such an expansion can be divided into three parts, the rst of which is that of the Galerkin method. The second part includes the in uence of the approximations of the right-hand side of the linear system, and the third one that of the approximations of the coe cients of the matrix (which require the evaluation of double integrals). Some additional notations about the spline space are called for. Let 0 be the characteristic function of the interval ( 1 2; 1 2) and m := 0 m 1 for all m 1. Then for m 1, m is a non-negative function with bounded support [ m+1 2 ; m+1 2 ], unit integral, of class Cm 1, piecewise polynomial of degree m, and so on. These and other properties can be reviewed in Aubin (1972) or Hsiao, Kopp and Wendland (1980). For i 2 Z, and h small enough, we de ne i to be the 1 periodic extension of the function m s zi h ; de ned in a unit length interval around zi. Obviously f igNi=1 is a basis of Vh. Let us nally consider a quadrature formula ZR m(t)f(t)dt ' Lm(f) := l X i= l cif(xi); where for all i xi = x i; ci = c i; and l is a certain positive or null integer. We will demand that the formula Lm be of degree m+ 1 at least, that means, that it be exact for polynomials of degree lesser than or equal to m + 1. However, the above demanded symmetry of nodes and coe cients implies that Lm will always have odd degree. For instance, for m = 0, we can take a single point formula (l = 0) L0(f) := f(0); which is exact in P1. For m = 1, the following three-point formula L1(f) := 1 12f( 1) + 5 6f(0) + 1 12f(1) is exact not only in P2, but also in P3. A three-point formula will also be su cient for m = 2, L2(f) := 1 8f( 1) + 3 4f(0) + 18f(1): 3 2 Approximation of the right-hand side The problem arising in the computation of the right-hand side is the evaluation of the integrals Z 1 0 u0(s) j(s)ds = h ZR u0(zj + th) m(t)dt; which, applying the quadrature formula, can be substituted by hLm(u0(zj + h : )) = h l X i= l ciu0(zj + hxi): This is equivalent to de ning an approximate linear form in Vh, (u0; rh)h := h N X j=1 rh;jLm(u0(zj + h : )); (3) where rh(s) = N X j=1 rh;j j(s): Because of their periodicity, the functions i are equal modulo N , so we can associate to rh 2 Vh the sequence frh;jgj2Z, with rh;j = rh;k for j k, modulo N . Consequently the summation in (3) can be done on any set of N consecutive indices. We denote [x] the largest integer which is lesser than or equal to a real number x. Proposition 1 For all u 2 C1, there exists a sequence of functions fungn 1 C1 such that for all rh 2 Vh, (u; rh)0 (u; rh)h = M X n=[m+1 2 ]+1 h2n(un; rh)0 +O(h2M+2)krhk 1 2 : Proof: For simplicity of notations we only prove the case m even. Since ZR m(t)dt = 1; we can write (u; rh)0 (u; rh)h = h N Xj=1 rh;j (zj; h) (4) where (s; h) := ZR m(t)0@u(s+ ht) l X i= l ciu(s+ hxi)1A dt: Let us remark that is in nitely often di erentiable in both variables, is 1 periodic in s, and satis es that for all s (s; h) = (s; h) (5) 4 because of the symmetry of nodes and coe cients. Moreover, since Lm is exact in Pm+1, we have that for k m+ 1 ZR m(t)0@tk l X j= l cixki1A dt = 0; from where @k @hk (x; 0) = 0; (6) for all k m+ 1. Denoting fk(s) := 1 (2k)! @2k @h2k (s; 0); by Taylor expansions in (4) and applying (5) and (6), it follows that (u; rh)0 (u; rh)h = M X k=m2 +1 h2k[fk; rh]h +O(h2M+2)krhk 1 2 ; (7) where [ ; rh]h := h N X j=1 rh;j (zj): (8) In Crouzeix and Sayas (1995) (see formula (29)) the existence of a sequence of real numbers fDjg, depending only on m (even or odd), such that for all 2 C1 and for all M [ ; rh]h = ( ; rh)0 + M Xl=1 h2lDl( (2l); rh)0 +O(h2M+2)krhk 1 2 (9) is proven. With uk(s) := k m2 1 X n=0 Dnf (2n) k n (s); the statement of the proposition is simply a consequence of (7) and (9). 3 Approximation of the bilinear form. Let us denote Jj := zj m+ 1 2 h; zj + m+ 1 2 h : Then the restriction to Vh of the bilinear form b can be written b(gh; rh) = N X i;j=1 bi;jgh;irh;j where bi;j := ZJi ZJj log jx( ) x(s)j2 i( ) j(s)d ds 5 and fgh;ig; frh;ig are the respective coe cients of gh and rh in the basis f ig of Vh. Our aim is to calculate approximately the coe cients bi;j. For simplicity of notations we will consider the usual sequence of coe cients frh;jgj2Z for each rh 2 Vh. From the periodicity of the functions i, bi;j can be de ned for all i; j 2 Z, modulo N . Let us de ne the index set I(N; j) := 8><>: fj p; j p + 1; : : : ; j + pg; if N = 2p + 1, fj p; j p + 1; : : : ; j + p 1g; if N = 2p: Thus, the bilinear form can also be written b(gh; rh) = N Xj=10@ X i2I(N;j)bi;jgh;i1A rh;j : With this choice of indices, ji jj N=2. Denoting F ( ; s) := log jx( ) x(s)j2 ( s)2 ; we can write, for all j and for all i 2 I(N; j), bi;j = ZJi ZJj F ( ; s) i( ) j(s)d ds+ h2(log h2 +Bi j;m); where Bk;m := ZR ZR log (k + (t u))2 m(t) m(u)dtdu; (10) or also, bi;j = h2 ZR ZR F (zi + ht; zj + hu) m(t) m(u)dtdu+ h2(log h2 +Bi j;m): (11) Let us notice that for all k, B k;m = Bk;m. Since F ( ; s) 2 C1(D), where D := f( ; s) : j sj 2 3g; we can easily apply the bidimensional quadrature formula proceeding from Lm and de ne an approximation of bi;j, for all j and for all i 2 I(N; j) i;j = h2 l X k;k0= l ckck0F (zi + hxk; zj + hxk0) + h2(log h2 +Bi j;m) (12) where the Bk;m are calculated exactly as shown in Hsiao, Kopp and Wendland (1980). That leads us to de ne an approximate bilinear form in Vh by bh(gh; rh) = N Xj=10@ X i2I(N;j) i;jgh;i1A rh;j: (13) 6 It is important to see that bh( : ; : ) is symmetric. If N is odd, then j 2 I(N; i) if and only if i 2 I(N; j), so the symmetry is obvious. If N = 2p, a more careful look is needed. If we de ne j+p;j by (12), then it is easy to see that j p;j = j+p;j, because of the properties of the coe cients and nodes of the quadrature formula Lm. Hence, we can de ne k := 8><>: 1 2 ; if k = p or k = p, 1; if jkj < p and rewrite (13), for the case N even, as bh(gh; rh) = N Xj=10@ j+p X i=j p i j i;jgh;i1A rh;j: (14) which proves that the bilinear form is symmetric. In Hsiao, Kopp and Wendland (1980) a way of de ning the coe cients i;j for all i; j (and obtaining the equality of coe cients when the indices are congruent modulo N) through the introduction of two discrete distances (i; j) := minfji j + kN j : k 2 Zg and (i; j; xn; xn0) := ji j + kN + xn xn0j; if (i; j) = ji j + kN j; is shown. Thus, i;j = h2 l X k;k0= l ckck0 log jx(zi + hxk) x(zj + hxk0)j2 h2 (i; j; xk; xk0)2 + h2(log h2 +B (i;j)): However, when N = 2p, (i p; i; xn; xn0) are not uniquely de ned. Nevertheless, again by the symmetry properties of the quadrature formula, both possible determinations of i p;i coincide, provided that we always choose the same determination of . For instance (i; j; xn; xn0) = ji j+kN+xn xn0j if (i; j) = ji j+kN j and (i; j) 6= ji j(k 1)N j. As in the previous section, we will restrict ourselves to the case m even for the sake of simplicity. At the end of the section the di erences will be pointed out for the case m odd. Lemma 2 If m is even or zero, there exists h1 > 0 such that for all h h1, jb(gh; rh) bh(gh; rh)j Chm+1kghk 1 2 krhk 1 2 ; for all gh; rh 2 Vh. Proof: Subtracting (11) from (12), we get that bh(gh; rh) b(gh; rh) = h2 N X j=10@ X i2I(N;j) gh;irh;j (zi; zj;h)1A (15) 7 where ( ; s;h) := ZR ZR(F ( +ht; s+hu) l X k;k0= lckck0F ( +hxk; s+hxk0) m(t) m(u)dtdu: The function belongs to C1(D [ h0; h0]) and satis es ( + 1; s+ 1;h) = ( ; s;h); ( ; s; h) = ( ; s;h); s+ 1 2 ; s;h = s 12 ; s;h and for all k m+ 1 @k @hk ( ; s; 0) 0: This last property is a consequence of the degree of the quadrature formula in the sense given in its de nition. Hence, jb(gh; rh) bh(gh; rh)j Chm+4 N X j=1 jrh;jj N X j=1 jgh;jj from where the result follows readily since h N X j=1 jrh;jj Ch 1 2krhk 1 2 for all rh 2 Vh. From the previous lemma, it follows that (at least for h small enough) the family of approximate bilinear forms bh satis es uniformly the Brezzi-Babu ska condition. Therefore (P h ) is uniquely solvable and its solution is uniformly bounded in norm H 12 by the data. The convergence of this method was proven in Hsiao, Kopp and Wendland (1980). We will deal with this aspect in the following section, before giving the asymptotic expansions of the error. Proposition 3 There exists ffngn m2 +1 C1 such that for h = 1=(2p + 1), and for all rh 2 Vh, b(Qhg; rh) bh(Qhg; rh) = M X n=m2 +1h2n(fn; rh)0 +O(h2M+2)krhk 12 : Proof: Let be the function de ned in the proof of the previous lemma and let k( ; s) := 1 (2k)! @2k @h2k ( ; s; 0); 8 for k m2 + 1, which satisfy k( + 1; s + 1) = k( ; s): Then, if faigi2ZN are the coe cients of Qhg in the usual basis of Vh, it follows bh(Qhg; rh) b(Qhg; rh) = h2 N Xj=10@ j+p X i=j p airh;j (zi; zj;h)1A = = M X k=m2 +1 h2k+20@ N Xj=1 j+p X i=j p airh;j k(zi; zj)1A +O(h2M+ 3 2 )krhk 1 2 ; because of the properties of . From properties of the operator Qh (see Sayas (1994)), we know also that there exists a sequence of numbers Cl independent of h and g such that ai = M Xl=0Clh2lg(2l)(zi) +O(h2M+2); uniformly in i 2 ZN . Thus, denoting for k m2 + 1 k( ; s) := k m2 1 Xl=0 Clg(2l)( ) k l( ; s); we have that bh(Qhg; rh) b(Qhg; rh) = = M X k=m2 +1h2k+20@ N Xj=1 j+p X i=j p k(zi; zj)rh;j1A +O(h2M+ 3 2 )krhk 1 2 ; (16) for all rh 2 Vh. Set now k;0(s) := Z s+ 12 s 12 k( ; s)d = Z 12 12 k( + s; s)d : Since for all ( ; s) 2 D, k( + 1; s + 1) = k( ; s); it follows that the functions k;0 are 1 periodic, and given that the functions k are smooth in D, so are k, from where k;0 2 C1. Remark now that h j+p X i=j p k(zi; zj) is the composite midpoint rule applied to compute k;0(zj). By the asymptotic expansion of the error in this formula h j+p X i=j p k(zi; zj) = k;0(zj) + M 0 X n=1 h2n k;n(zj) +O(h2M 0+2); (17) 9 where k;n(s) := 2n Z 1 2 12 @2n @ 2n k( + s; s)d and 2n are some real coe cients. Note that k;n are also smooth and 1-periodic. Applying (17) with M 0 = M k to (16), and denoting 'n(s) := n m2 1 X k=0 n;n k(s); it follows that for all rh 2 Vh, bh(Qhg; rh) b(Qhg; rh) = M X n=m2 +1h2n['n; rh]h +O(h2M+2)krhk 1 2 ; where the discrete product [ : ; : ]h is de ned by (8). Finally, (9) gives the result. Proposition 4 There exists f ngn m2 +1 C1 such that for h = 1=2p, and for all rh 2 Vh, b(Qhg; rh) bh(Qhg; rh) = M X n=m2 +1 h2n( n; rh)0 +O(h2M+2)krhk 1 2 : Proof: Let k be the same functions as those de ned in the proof of Proposition 3, and let us remark that since k s+ 1 2 ; s = k s 12 ; s ; then k s+ 12 ; s = k s 12 ; s for all s and consequently k(zj p; zj) = k(zj+p; zj); for all j. Hence h j+p 1 X i=j p k(zi; zj) is the composite trapezoidal rule applied to compute k;0(zj) (being this function the same as in the previous proposition). The remainder of the proof of Proposition 3 is valid, changing the coe cients appearing in the asymptotic expansion of the error for the midpoint rule by those for the trapezoidal rule. Remark. For m odd, we have that for all rh; gh 2 Vh jb(gh; rh) bh(gh; rh)j Chm+2kghk 1 2 krhk 1 2 ; and consequently existence and uniqueness of solution follows. Propositions 3 and 4 also hold here, changing only the lower limit of the sum, which now is n = (m + 3)=2. The expansions for N even and odd are also di erent. 10 4 Asymptotic expansion of the error In Hsiao, Kopp and Wendland (1980), the following convergence result is proven for m = 0; 1 and 2, although it can be extended to arbitrary m. Theorem (convergence). Let u0 2 Hm+2, g be the unique solution of (P ) and g h the unique solution of (P h ). Thenkg h gks Chm+1 skgkm+1 (18) for all s 2 [ 1 2;m], if m 1, and for s = 0 if m = 0. Let us simply remark that the original result in that paper is stated for a slightly modi ed problem. Propositions 1, 3, 4 and the decomposition of bh(Qhg g h) given in (2) prove the following result. Proposition 5 For all positive integer M , there exist functions ek, as smooth as desired, such that bh(Qhg g h; rh) = M X k=2[m+2 2 ]hk(ek; rh)0 +O(hM+1)krhk 1 2 (19) for all rh 2 Vh, if h = 1=(2p + 1). A similar expansion also holds when h = 1=2p. Let us consider the Sobolev space Wm;1 := fu 2 L1 : u(j) 2 L1; 1 j mg; with the essential supremum norm k : km;1. Theorem 6 For all positive integer M , there exist functions fk, as smooth as desired, such that kg h Qhg M X k=2[m+2 2 ]hkQhfkkm;1 ChM+1 (20) if h = 1=(2p + 1). A similar expansion also holds when h = 1=2p. Proof: The proof of Theorem 2 of Crouzeix and Sayas (1995) is valid here, applying now Proposition 5 and the properties of the family of approximate bilinear forms bh. Since the convergence theorems give the same orders of convergence in the same `strong' norms for the Galerkin and the Galerkin collocation methods, the other results of this section do not di er essentially from those given for the Galerkin method. Corollary 7 Let h be the set of the midpoints of the grid. Then for all 0 l m, there exists a sequence of constants fFk;l(z) : z 2 hg, independent of h = 1=(2p + 1), such that sup z2 h j(g h)(l)(z) g(l)(z) M X k=2[m l+2 2 ]hkFk;l(z)j = O(hM+1) 11 for all positive integer M . A similar expansion also holds for all h = 1=(2p).Let h be the set of the nodes of the grid. Then for all 0 l m 1, there exists asequence of constants fFk;l(z) : z 2hg, independent of h = 1=(2p + 1), such thatsupz2 h j(gh)(l)(z) g(l)(z)MXk=2[m l+22 ]hkF k;l(z)j = O(hM+1)for all positive integer M . A similar expansion also holds for all h = 1=(2p).Finally, what is said about superconvergence in L1 norm for the Galerkin method istrue here too. Likewise it is possible to build an interpolate solution of higher order whenm is even as was pointed out there.5 Numerical resultsWe have tested the fully discrete Galerkin collocation method with piecewise constantfunctions (that is, m = 0), with the midpoint rule as numerical quadrature (see examplesgiven at the beginning), on some problems whose exact solution is known. We show heresimply an example where it can be seen thatkQhg ghk1 Ch2;(21)and that Richardson extrapolation can be applied as a means of accelerating the con-vergence. Remark that from the last corollary it follows that we have `true' asymptoticexpansions (with all coe cients independent of the numerical parameters) when restrict-ing ourselves to a nite set of points, all belonging as midpoints to all meshes underconsideration. The number of points where extrapolation is to be applied restricts thesize of the di erent grids where the solutions are calculated. As the simplest test we showhere results on a single point.Let be the ellipse (x=R)2 + y2 = 1. Then with datau0(x; y) := x3 3xy2; (x; y) 2 ;the solution of the variational problem is the product ofq(x; y) = 3(R + 1)22R2xq x2R4 + y2 (R 2(R + 1)y2);where x = (x; y) is substituted by the usual parametrization of , by the norm of theparametrization jx0(s)j2. We choose R to be equal to 2:Figure 1 shows the midpoint convergence of the solution computed for N = 10; 20;30; : : : ; 400. We have plotted there the decimal logarithm of the number of intervalsagainst log10(EN ) where EN stands for the maximum error in the midpoints, that is,(21). Applying linear regression to that set of points we obtain an estimate of the orderof the method. The slope of the regression line is in this case 1.962.12 Figure 1: convergenceWe next show a table of the errors obtained by extrapolation when calculating thesolution in (R; 0). The solution has been computed for N = 20; 40; 80; 160 and 320; that isin a Romberg set, making midpoints coincide. Richardson extrapolation has been appliedto solutions in consecutive grids. The table shows the errors: columns show di erentstages of the extrapolation process, rows show errors in correspondence with the nestgrid used for each computation.N = 20 1:35N = 40 3:53 10 1 2:0 10 2N = 80 9:04 10 2 2:65 10 3 1:73 10 4N = 160 2:28 10 2 3:66 10 4 3:93 10 5 3:04 10 5N = 320 5:75 10 3 4:86 10 5 3:31 10 6 9:09 10 7 4:44 10 8Numerical results show how Richardson extrapolation gives a remarkable increase ofprecision, in addition to faster convergence. Other numerical tests for this method can beseen in Sayas (1994).We wish nally to indicate that when solving some di erential problems by a bound-ary element method, the approximate solution calculated by some postprocessing of thenumerical solution of the BEM can inherit the existence of asymptotic expansions andthus the possibility of applying directly extrapolation.ReferencesAubin, J.P. (1972): Approximation of elliptic boundary-value problems, Wiley-Inter-science.13 Crouzeix, M. and F.J.Sayas (1995): Asymptotic expansions of the error of BoundaryElement Methods I: spline Galerkin methods, submitted for publication.Hsiao, G.C, P.Kopp and W.L.Wendland (1980): A Galerkin collocation Method for someIntegral Equations of the First Kind, Computing 25, 89-130.Hsiao, G.C, P.Kopp and W.L.Wendland (1985): Some applications of a Galerkin-Col-location Method for Boundary Integral Equations of the First Kind, Math. Meth.in the Appl. Sci. 6, 280-325.Le Roux, M.-N. (1974): Resolution numerique du probleme du potentiel dans le planpar une methode variationelle d'elements nis, These de troisieme cycle, Universitede Rennes, France.Sayas, F.J. (1994): Asymptotic expansion of the error of some boundary element meth-ods, Tesis doctoral, Universidad de Zaragoza, Spain.Sloan, I.H. (1992): Error analysis of boundary integral equations, in Acta Numerica,Cambridge University Press.14